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ABSTRACT 

A method of deriving and using merging history trees of dark matter galaxy haloes 
directly from pure gravity N-body simulations is presented. This combines the full non- 
linearity of N-body simulations with the flexibility of the semi-analytical approach. 
Merging history trees derived from powe r-law initial pe rturbation spectrum sim- 



ulations (for indices n = —2 and n = 0) by Warren et al. (1992) are shown. As an 



1 INTRODUCTION 



As q uantitativel y pioneered by th ose such as Hoyle (1953) 



Silk (1977) and Rees & Ostrikei (1978) and more recently 
noted by authors such as Frenk etal. (1995), hierarchical 
galaxy formation models in f2o = 1 cold dark matter (CDM) 
universes typically combine assumptions on up to six dis- 
tinct physical processes: (1) the non-linear growth phase of 
matter density peaks (known as "haloes"), (2) cooling gas 
dynamics, (3) star formation, (4) star-to-gas energy feed- 
back, (5) stellar evolution, (6) galaxy mergers. In principle, 
if there are more free parameters describing these processes 
than i ndependent observational galaxy statistics, then the 



example of a galaxy formation model, these are combined with evolutionary stellar 
population synthesis, via simple scaling laws for star formation rates, showing that if 
most star formation occurs during merger-induced bursts, then a nearly flat faint-end 
slope of the galaxy luminosity function may be obtained in certain cases. 

Interesting properties of hierarchical halo formation are noted: (1) In a given 
model, merger rates may vary widely between individual haloes, and typically 
20%~30% of a halo's mass may be due to infall of uncoUapsed material. (2) Small 
mass haloes continue to form at recent times: as exp ected, the existence of young, low 
redshift, low metallicity galaxies (e.g., [zotov ct al. 1997) is consistent with hierarchi- 
cal galaxy formation models. (3) For n ~ —2, the halo spatial correlation function can 
have a very high initial bias due to the high power on large scales. 

Key words: Methods: numerical - Galaxies: formation - Cosmology: theory - Galax- 
ies: luminosity function, mass function - Galaxies: interactions - Galaxies: stellar con- 
tent. 



lation synthesis for process (5). Since each of these models 
have problems explaining at least some of the observations 
means, the models are better constrained than might have 
been hoped for. 



observations should provide little constraint on galaxy for- 
mation "recipes" . Fortunately, the contrary is presently the 
case for the "semi-analytical ab initio" models which make 
various analytical estimates of process (1), combine semi- 
empirical and simple scaling parametrisations to represent 
processes (2)- (4) and (6) and use evolutionary stellar popu- 



These models can be considered to be semi-analytical 
because rather than calculating what is possibly the most 
important process, the non-linear formation and merging 
history of collapsed objects [process (1)], via N-body simula- 
tions, various statistical analytical approximations are used. 
The model s of Lacey et a l. (1993) use an approximation de- 
veloped in j^^accy fc Sills (1991) from the BBKS peaks for- 
malism (Bardccn ct al. 1986), Kauffmann et al. (Kauffmann 



& White 1993 



Kauffmann et al 



1993) use a probabilistic 
method ( Bower 1991) based on the Press-Schecter formal- 
ism (e.g.. Press fc Schechtet| 1974; White fc Frenk 1991 ) and 



excur sion set mas s function calculations (Bond et al 
while 



Cole et al 



1991) 

(1 994) use a spat ially quantised "block" 
method described in Cole fc Kaiser (1988). 
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Further semi-analytical developments include those 
adding spatial auto-co rrelation inform ation to a Press- 
Schechter form alism (Mo fc White 1996 ) or to the 
"block' 



^ ^ 3}- 

model (Rodrigues fc Thomai^ 1996; [Nagashima &] 

Gouda[]l997) and a technique of separately treating global, 

weakly non-linear and local, strongly non-lin ear dynamics 

(the "peak-patch" formalism. Bond & Myers 1996). 



Each of the models which has been compared to obser- 
vational statistics has difficulty in simultaneously explain- 
ing the flatness of the present-day (surface-brightness lim- 
ited) galaxy luminosity function (e.g., Loveday et al. 1992), 
the steepnes s of th e faint galaxy counts (e.g., Tyson & 
Seitzei^988; Tyson 1988), the shape of the moderately faint 
galaxy (j? ^ 23) s pectroscopic redshift distributions (e.g. 



CoUess et al] 1990; poUess et al] 1993), the TuUy-Fisher rela- 
tion and the colour distributions of prese nt-day galax ies, in 
a CDM Qo = 1 universe. Even though the |Cole et al.| (1994) 



model is better than the previous models in allowing big 
z = galaxies to be at least as red as higher z galaxies, it 
shares the problem of the other r nodels in lacking b ig red 

model 



Kauffmann et al 



ellipticals. It also shares with the 
the problem that if the large number of small haloes pre- 
dicted by CD M models at z ^ foUow the IR TuUy-Fisher 
relation (e.g.. Pierce & TuUy 1992), then the slope of the 



faint end of the general galaxy luminosity function should 
be steeper than that estimated locally (e.g., Loved ay etal. 
1992). Chan ging the cosmology in the Cole et al. models 
(Heyl et al. 1995: low Ho, low Qq, non-zero Ao and CHDM 
models) is insufficient to match the observations. Another 
way of allowing these models to fit the observations is to 
make a strong assumption for process (6) — to suppose that 
galaxies can merge as fast as galaxy haloes merge, or even 
faster — but si mple presen t-day constraints on the products 
of the mergers (Dalcanton 1993) and the relative w eakness of 
the faint galaxy angular auto-correlation function ( Roukema 
& Yosh|^ 1993) strongly restrict this possibility. 

In order to avoid problems which may be due to the 
approximation of non-linear gravitational collapse and evo- 
lution by the semi-analytical techniques mentioned above, 
an alternative technique is to calculate both processes (1) 
and (2) from first principles in numerical N-body simula- 
tions, folding in the other physical processes as simple scal- 
ing formulae or using stellar popul ation synthesis fo r (5). 
Several authors (e.g., Evrard 1988; [Nav arro fc Benz 1991; 



Cen & Ostrikei 1992; Umemura et al. 1993; Steinmetz & 
Mullei[]l995) have experimented with these techniques, but 
resolution limits on present-day computers mak e the results 
hard to interpret. For example, Weinberg et al. (1997) point 



out that although low resolution gravito-hydrodynamic sim- 
ulations suggest that a UV photoionisation background can 
suppress galaxy formation (by heating the gas so that it is 
unable to cool and form stars) , higher resolution simulations 
show that this is a numerical artefact: the higher resolution 
simulations show little sensitivity to either the details of 
photoionisation or star formation. 

In this article, rather than claiming a global "recipe" 
for galaxy formation, our primary purpose is to concentrate 
on process (1) in a way complementary to that of other 
techniques. This is unlikely to be sufficient to solve all the 
observational conflicts. On the contrary, this method should 
increase the ability of modellers to verify the extent to which 



model predictions are sensitive to the precision of modelling 
of gravity. 

The method presented here is to derive merging history 
trees of dark matter haloes directly from N-body simula- 
tions. Rather than just investigating virialised haloes for a 
particular dark matter model (e.g., CDM), (a) both n = —2 
and n = initial perturbation spectrum simulations (where 
n is the index of the power spectrum) are examined, and (b) 
since the halo-to-galaxy relationship may be more complex 
than a simple one-to-one mapping, two significantly different 
density thresholds are used for halo detection. This reveals 
the sensitivity of halo merger history trees and halo statis- 
tics to these param eters. The N-body simulations used are 
presented in 



in 52.1.2 



2.1.1 



the choice of a group-finding algorithm 
and the defining criterion and a lgorithm for calcu- 



lating the merging history trees in §2.1.3 



Properties of the haloes detected are discussed in §2.2 



In particular, the resulting merging history trees are pre- 
sented in graphical form in §2.2.3, enabling patterns of halo 
merging calculated from fully non-linear simulations to be 
visualised directly. 

If processes (2)-(5) are simple enough, and if process 
(6), galaxy merging, corresponds in a one-to-one way with 
halo merging, then these halo merger history trees would 
lead directly to galaxy merger history trees. We therefore 
examine an example application of the merger history trees 
by making minimal assumptions for processes (2)- (4), using 
stellar evolutionary population synthesis for process (5) , and 
for process (6), assuming maximal galaxy merging (every 
halo merger corresponds to a galaxy merger) . §p.l| presents 
(6) + (3) merger-induced star formation and j |3.2| explains 
how process (5) is modelled. 

In order for these processes to have an effect on the 
luminosity function, an option is considered in which each 
merger induces a burst of star formation, scaled according 
to the appropriate halo and gas masses and the dynamical 
time scale. Apart from this star formation rate option, we 
do not explore parameter space for non-gravitational pro- 
cesses in this paper; we merely adopt simple observationally 
normalised scaling laws. Resulting luminosity functions are 
presented in §3.3. 

Applications of N-body derived halo merger trees with 
more complex assumptions for processes (2)- (6) are of course 
possible, and indeed to be welcomed. The galaxy formation 
"recipe" explored here is only one simple example. 

Cosmological conventions adopted for this paper are a 
Hubble constant of Ho = 50kms~^Mpc~^, comoving units 
(at t = to) and an fio ~ TO, A — universe is assumed, 
except where otherwise specified. 



2 HALO MERGER HISTORIES (GRAVITY) 
2.1 Method 

2.1.1 N-body Models of Matter Density 

The non-linear gravitational evolution of matter density is 



modelled by N-body cosmological simulations run by War 
(1992). These simulations use a 128^ initial particle 



ren et al. 



mesh, of side length 10 Mpc. [The simulations analysed here 
are for power law initial perturbation spectra (n — —2 and 
n = 0), so this is simply a default choice for the scaling of 
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units. This default scaling is used hereafter except where 
otherwise specified.] Particles are placed on this mesh, mak- 
ing a cube of ~ 2x10^ particles. 

An initial perturbation spectrum is imposed on this 
cube by Fourier transforming the initial complex amplitudes 
from the perturbation spectrum and using the Zel'dovich 
growing mode method (Warren et al. 1992) on this Fourier 



transform and the 128 particle mesh. The amplitude of the 
perturbation spectrum is chosen such that linear perturba- 
tion growth implies that {5M/M)(r = 0-5ft"^Mpc) = 2-0 
at z = 0, where {5M /M){r) is the r.m.s. value of the excess 
m ass (ov er uniform density) in spheres of radius r ( Warren et 
al 



199; 



This choice ensures that the haloes which collapse 
are about the same size for different values of n, so that the 
dependence on n of properties of halo dynamics — or merging 
histories — can easily be studied. The absolute normalisation 
of the spatial correlation function of the haloes cannot be di- 
rectly interpreted in terms of observational quantities. The 
relative amount of power on different scales (or slope of the 
spatial correlation function), and the halo detection thresh- 
old, are the parameters which may affect the rates and ways 
in which haloes merge with one another. 

The initial cube of perturbed particles is trimmed to a 
sphere, i.e., particles more than 5000 kpc from the centre of 
the cube are removed, resulting in a sphere of ~ 1-lxlO'' 
particles. 

This is the n evolved forwa rd gravitationally via a tree- 
code (e.g., see Barnes & Hut 1986), initially with roughly 



logarithmic time steps up to f = 0-3 Gyr, after which equal 
time steps of 0-03 Gyr are used. Every hundredth time step 
is stored on disk; these are the time steps available for halo 
analysis (hereafter "time stages" ) . A vacuum boundary con- 
dition is used and the softening parameter is 5 kpc (proper 
units). 

2.1.2 Group-Finding Algorithm 

The simulation data are searched for density peaks at each 
time step by an algorithm which uses the "oct-tree" method 
to find all overdense regions without overwhelming computer 
memory, followed by an iterative means of joining together 
contiguous overdense regions. 

Alternative group-finding methods which could be 
used include the "friends-of-friends" (FOF) algorithm 
(e.g.. White et al 



ren et al. 



1987), the algorithm used by War- 



(1992) or the DENMAX algorithm (Gelb & 



Bertschingei 1994). 

The FOF group-finder has the advantage of low mem- 
ory requirements and an obvious relation between the mean 
particle separation and the group-finding resolution, but has 
the disadvantages that if the link parameter I is too low, then 
low density haloes — or the low density envelopes of haloes — 
are missed, while if / is higher, small but distinct haloes may 
be erroneously join ed together as single objects. 

Warren et al.'s (1992) method, based on the accelera- 
tions of individual particles, and the DENMAX algorithm, 
which includes a de-binding procedure to separate haloes 
which are only temporarily close to one another, are both 
more physically motivated than FOF. However, for a first 
investigation of the use of N-body generated merging his- 
tory trees in galaxy formation models, the use of the simple 
method outlined below seems prudent. Since two different 



density detection thresholds are used, the implications of 
having either a low or a high fixed density threshold (which 
are similar to the cases of high or low I respectively in FOF) 
can be seen. For further development, it would certainly be 
useful to consider use of a more complex algorithm such as 
DENMAX. 

Details of the method are as follows. 

Conceptually, a cube concentric to the sphere of par- 
ticles, having as side length the diameter of the sphere, is 
divided into eight equally sized subcubes. Any of these sub- 
cubes containing more than one particle is itself subdivided 
into eight subcubes. By not subdividing cubes with only 
one or zero particles, computer memory is not wasted on 
analysing "empty" space. The subdividing process is iter- 
ated to a depth of nigygjg levels below the original cube, 
unless at some level all the cubes have one or zero particles 
in them, in which case subdividing stops (this would happen 
at njg-^rgjg = 8 for this l-lxlO^'-particle model for a uniform 
particle distribution). The side length of the smallest cube 
is 174 kpc and 20 kpc for wjgygjg = 6 and n[g.^,gjg — 9 respec- 
tively at 2 = 0. 

The "primary" list of density peaks is then simply the 
list of each cube at the deepest level (i.e., of size 2""''=^"''' 
times the simulation sphere diameter) whose density is at or 
above r^hresh times the mean density. The list of particles 
in each of these peaks is recorded. 

The results presented here are for r^hresh ~ 5 and 
''thresh = 1000. For a flat rotation curve of the Galaxy 
of 220kms~^, the cumulative mass to a radius r is Af(< 
r) (X r, and the density is p{r) — 1-2x10^ pcr~^ for Ho as 
above and r in kpc. So, detection densities of r^-j^j-ggj^ — 5 
and f^hresh ~ 1000 correspond to the total (baryonic plus 
nonbaryonic) matter density at galactocentric distances of 
about 1500 kpc and 110 kpc respectively. The latter is a rea- 
sonable value for the halo radius, but the former is several 
times greater than the largest radii claimed for the halo of 
the Galaxy. 

The key to a simple way to join together contiguous 
"primary" peaks is to order the primary density peak list 
by mass, from largest to smallest, so that each "secondary" 
peak can be created by starting from its nucleus (densest 
region) and successively joining on regions of lower and lower 
density which are adjacent to the region which has already 
been aggregated.^ 

The parameters used to decide on adjacency are the ra- 
dius r (from centre to outermost particle) of each secondary 
peak and an "incremental radius", rj^^, defined as IT times 
half of the largest diagonal of the small cube used in finding 
the primary density peaks. A primary peak is considered ad- 
jacent to a secondary peak if its centre is within r -I- rjj^g of 
the secondary peak. The radius r is re-evaluated each time 
a primary peak is joined to a secondary peak. The nucleus 
of the first secondary peak is the first (i.e., most massive) 
primary peak; each following secondary peak starts with the 
most massive primary peak not previously included in a sec- 
ondary peak. 

The final secondary peak list, corresponding to all sepa- 
rated regions inside isodensity contours of ^'thresh times the 



* The idea to order the "primar y" peaks bv rn ass was inspired 
by the group-finding algorithm of Warren et al. (1992). 
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mean density, is hereafter simply termed a peak list, since 
the list of primary peaks is not of astrophysical interest. The 
members of this list are considered to be dark matter haloes. 



2.1.3 Creation of History Tree 

A peak (halo) merging history tree is obtained as follows. 

Peak lists for a series of time stages of the model are ob- 
tained by the algorithm just described, each obtained with 
the same values of ni^^^i^ and '"thresh- For each pair of suc- 
cessive output times, ti,ti+i, the peaks at the two times are 
compared. Two arrays, ai, Oi+i, each with as many elements 
as the number of particles in the simulation, are created. For 
each element j of array at, the integer k identifying the peak 
that particle j is a member of is assigned to ai{j), where this 
is a peak according to the peak list for ti. The array a^+i 
is evaluated in the same way using the peak list for ti+i. 
If the particle is not a member of any peak, a null value is 
assigned. A simultaneous sort is performed on arrays ai and 
ai+i, permuting both in the same way in order that a^+i is 
a non-decreasing arithmetical sequence. 

The result is that with the new ordering of ai and 
ai+i, (1) a peak k at ti+i is represented by a contiguous 
list ai+i{j) to ai+i{j') (each containing the peak number k) 
and (2) ai{j) to ai{j') (for the same j, j') represent the same 
particles and contain values (fei, fe2, ...) indicating the peaks 
(at ti) of which the particles were members. In other words, 
the peak membership at ti of particles in a single peak at 
ti+i is listed in ai{j) to ai{j'). 

For any peak at ti+i, if more than 50% of the particles 
in any of the peaks at ti are present in the peak at ti+i, then 
the peak at ti+i is considered a "descendant" of the peak at 
ti and the peak at ti is a "progenitor" of the peak at ti+i. 
These links are represented by appropriate arrays. Due to 
the nature of this algorithm, no peak can have more than 
one descendant, though it can certainly have more than one 
progenitor, which is allowed for by using what are, in effect, 
pointers. 

By applying this comparison across each pair of suc- 
cessive times ti,ti+i, a representation of the peak merging 
history is obtained. 

2.2 Results 

The method described above has been applied to both an 
n = and an n = —2 power law initial perturb ation spec- 
trum N-body model (labelled "nOb", "n-2b" by Warren et 
1992). Table |l| shows redshifts and cosmological times 
for the output timesteps for these two models. The negative 
redshifts correspond to future times according to the default 
time scaling. If the time unit chosen were different to the de- 
fault, then these latter time stages could be moved into the 
past or the present. 

Haloes are detected in both of these models at thresh- 
olds of both r^]-|j.ggjj = 5 and ^'thresh ~ 1000 times the mean 
universe density (for anf2o = l-0,/i = 0-5,Ao = universe.) 
The former density threshold for detection should result in 
haloes which have only just turned around from following 
the smooth Hubble flow, while the latter should result in 
well-virialised haloes. These thresholds should span most 
cases of interest. Statistics of the haloes detected and com- 
parison of these to results from other methods are presented 



Table 1 . Parameters of Time Stages Used 



redshi ft 


*(tjyr) 


timestep 






n = 


n = -2 


11-2 


0-31 


40 


(15) 


3-2 


1-51 




55 


1-5 


3-3 


140 


115 


0-62 


6-3 


240 


215 


0-25 


9-3 


340 


315 


039 


12-3 


440 


415 


-0-10 


15-3 


540 


515 


-0.203 


18-3 


640 


615 



Table 2. Number of haloes found for the different power spectra 
and detection thresholds. 



t(Gyr) 



''thresh 



^thresh 



1000 



: n : 



: n -- 



0- 3 3959 - 238 

1- 5 - 2086 - 412 
3-3 1539 1890 4214 1421 
6-3 1053 1121 2695 1516 
9-3 836 804 2121 1176 
12-3 712 637 1891 923 
15-3 629 492 1674 790 
18-3 590 433 1524 672 



2.2.1 and 



displayed in §2.2.3 



2.2.2 ; examples of halo merging histories are 
the spatial two-point autocorrelation 



functions of the haloes are discussed in §2.2.4; and a sug- 
gestion for further development by inte rpolat ion of merger 
times between time steps is outlined in ^2.2.5. 



2.2.1 Basic Halo Statistics 

The numbers of haloes are shown in Table ^. In the n = 
—2 model, matter has not yet collapsed into haloes ai t — 
0-31 Gyr, so the t — 1-5 Gyr time step was used instead. 

The reality of these haloes is verified visually by recti- 
linear projections of a sample of the points for each halo, 
by radial particle count profiles and by an interactive pro- 
gram which plots a sampling of all the points on a computer 
screen. The program offers optional colouring of a range of 
haloes in a colour different to both the particle colour and 
the black background and allows real time rotation of the im- 
age in order to give an intuitive feel for the three-dimensional 
shape of the data. Examples of haloes are given in Fig. |l|, 
which shows the radial particle count profiles for four of the 
biggest haloes (by number) detected using r^j^j-ggj^ = 5 in the 
final time stage of the n = — 2 model. The profiles are simply 
numbers of particles per spherical shell, so the rapid decrease 
to zero shows that the density falls off faster than r^ . Note 
that one profile has two prominent maxima, neither at r = 0. 
This is because, as a closer visual investigation of the haloes 
shows, a small proportion of the "haloes" are in fact fairly 
close binary haloes rather than single haloes. These binaries 
are usually quite uneven in size, so consideration of the halo 
as a single halo is still a good approximation. 

As noted above, use of a group-finder such as DEN- 
MAX (Gelb fc Bertschinger 1994) would be an elegant way 
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Figure 1. Halo profiles: numbers of particles in spherical shells 
for some of the most massive haloes ( "peak" number = 1 — 3, 12 as 
labelled) of time stage 615 in the n = —2 model for »"^iji.ggjj = 5. 
The total number of particles in each peak is labelled "ninpk" . 



Table 3. Statistics of fraction of halo at time stage listed here 
contained in halo at following time stage. 



Figure 2. Mass functions for the n = —2 model for ''thresh ~ ^' 
Thick lines are for this work at th e time stages labelled. For com- 
parison with ot her results (§2.2.2), the Press & Schechter formula 
(Laccy fc Cole 1993: thin lines; 5^0 = ^vir = 1-686) and the mod- 



els of Laccy fc Silk (1991) (thin solid line) and of Bond fc My- 



(1996) (hollow and solid circles for spherical and ellipsoidal 
internal dynamics resp.) for a CDM initial fluctuation spectrum 
are shown. 



t(Gyr) 


''thresh — ^ 


'"thresh 


= 1000 


n = 


n = -2 


n = 


n = -2 


0-3 


96 ± 8% 




90 ± 11 % 




1-5 




87 ± 15% 




76 ± 12 % 


3-3 


92 ± 9% 


84 ± 14% 


78 ± 12 % 


74 ± 12 % 


6-3 


92 ± 9% 


81 ± 14% 


79 ± 12 % 


71 ± 11% 


9-3 


93 ± 8% 


82 ± 13% 


81 ± 12 % 


72 ± 10 % 


12-3 


93 ± 8% 


82 ± 13% 


82 ± 11 % 


72 ±11% 


15-3 


94 ± 8% 


84 ± 12% 


84 ± 11 % 


74 ± 10 % 


18-3 











to avoid the adoption of "binary" haloes as single haloes. Al- 
ternatively, a positional proximity based group-finder could 
separate bound overlapping bound groups of particles by 
requiring an additional layer of analysis, such as detection 
of "haloes within haloes" , in which either f^-j^j.^gj^ or I (for a 
friends-of-friends group-finder) would have to again be spec- 
ified. 

As described in 



.1.3 



for each time stage, a halo at 
time ti is considered to merge into a halo at the following 
time stage ti+i (or retain its identity) if and only if more 
than 50% of the particles of the halo at time ti axe present in 
the halo at time ti+\. Table ^ shows the means and standard 
deviations of the fraction of a halo at time ti present in a halo 
at time ti+i. By definition, these fractions are constrained 
to be greater than 50%. 

Figs ^-|^ show the mass functions of these haloes and 
comparable mass functions from various semi-analytical 
methods. Our N-body derived mass functions are inter- 
preted in terms of merging rates in the following para- 
gra phs a nd compared with other mass function calculations 
in 



2.2 




10'" 10" 10" 

mass of peak m 



The mass function figures suggest that with the detec- 
tion threshold of r^hresh ~ 5, for either spectral index the 



Figure 3. Mass functions for n = model for '"thresh 
Curve styles are as for Fig. ^. 



overall halo merging rate from t = 3-3 Gyr (i.e., for z ~ 1-5), 
to the present is little more than about 3 — 10 for galaxies 
below about 10^'' M0. While this merging goes on, the num- 
ber of large haloes in the largest bin in the n — model 
increases somewhat until the last time step. Depending on 
the average number of small haloes which merge into a single 
large one, the increase in the number of large haloes would 
appear at first sight to be explained by the decrease in the 
number of smaller ones, consistent with a merging ratio of 
about 3— 10. Though the high mass end of the n = —2 mass 
function is noisy, similar interpretation could be made. 

These plots show a significant dependence on detection 
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Figure 4. Mass functions for n = —2 model for J'thi-gsh ~ 1000. 
Curves are as for Fig. ^ except that <5cO = '^^yiv formula 
for the 1-5 Gyr and 3-3 Gyr time steps. 



Table 4. Fraction of haloes which have no descendants at fol- 
lowing time stage. 



*(Gyr) ''thresh = 5 ^thresh = ^O"" 
n = n = —2 n = n = —2 



0-3 


5% 




11% 




1-5 




15% 




32% 


3-3 


8% 


24% 


32% 


40% 


6-3 


8% 


30% 


26% 


49% 


9-3 


7% 


29% 


23% 


46% 


12-3 


7% 


24% 


20% 


44% 


15-3 


4% 


21% 


15% 


36% 


18-3 











Table 5. Numbers of original haloes which end up in a halo 
detected at the final time stage (meanist. deviation). 



'•thresh " = " n = ^2 

5 7-4 ± 20-7 5 ± 16-9 
1000 3-2 ± 6-5 2-6 ± 6-2 




JQIU jgll 

mass of peak m Mg 



Figure 5. Mass functions for n = model for ''thresh ~ 1000 
Curves are as for Fig. ^ except that S^o = 
for the 1-5 Gyr time step and <5cO = 25- 
step. 



5<5.yjj, in the PS formula 
for the 3-3 Gyr time 



threshold and a weak dependence on n. For the haloes de- 
tected at r^jjj.ggjj — 1000, the merging is much weaker than 
for r^^hresh ~ ^- In a given simulation, objects detected at 
the higher threshold consist of the dense cores of the objects 
detected at the lower threshold. Hence, a simple explanation 
for the weaker merging is that if the low density envelopes 
merge, the cores don't necessarily do so, but if the cores 
merge, the large low density envelopes are almost certainly 
going to merge. 

However, some simple statistics show that the merging 
history is not as simple as an 7V:1 ratio applying equally to 
all haloes. 

Table ^ shows the fraction of the haloes at each time 
stage that have no descendants, i.e., the fraction of the 



haloes for which no more than 50% of their particles appear 
in any single halo at the following time stage. The fact that 
these are nonzero (from about 5% for n — 0, r^jjj-ggjj = 5 to 
30% — 50% for n = —2, r^hi-gsh = 1000) shows that many 
haloes are destroyed in the sense that more than 50% of their 
particles may have been pulled into an "atmosphere" of a 
large halo at a density lower than the threshold density or 
possibly thrown out of the halo or pulled into another halo. 
This means that the halo number density does not only de- 
crease by merging, it also decreases by this halo destruction. 
For example, if the overall number ratio between two time 
stages is 4:1, but one in four haloes terminates, then the 
underlying ratio of haloes actually merging is only 3:1. Of 
course, this distinction is dependent on the definition of halo 
detection identity as described above. 

More direct statistics are those of the histories of the 
haloes detected at the final time stage. The mean (and stan- 
dard deviation) of the overall number of haloes which col- 
lapse to above the threshold density (either at the first time 
stage or at a later time stage) and end up in a final halo is 
shown in Table ^ 

While these mean values are in the range 3 — 10 esti- 
mated above, the standard deviations show that many final 
haloes come from as many as 20 or more original haloes. 
In fact, the maximum number of haloes that any final halo 
originates from is 233 for the n — model and 259 for the 
n = —2 model (for r^hresh = 5). For r^j^j-ggj^ — 1000, the 
overall rate is lower, and the maximum numbers of haloes 
per any final halo are 88 and 95 for n — and n — —2 
respectively. 

Table 6. Mass fraction of a final halo which comes from matter 
directly accreted from the "background" (meanist. deviation). 



''thresh 


n = 


n = -2 


5 


32 ± 26 


23 ± 28 


1000 


36 ± 25 


29 ± 32 
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As already suggested by the number of haloes which ter- 
minate, throwing matter back out into the background, the 
amount of matter which "rains" or accretes onto haloes di- 
rectly r ather than first collapsing into smaller density haloes 



The time evolution of the PS functions, both in nor- 
malisation and shape, is similar to that of our haloes. This 
is best seen in Fig. ^, where the decrease in normalisation 
in pr oportion to (1 -I- z) [expected from Eq.(2.1) of Lacey & 



is non-i .egligible. Typically 30% of a final halo's mass comes 



from matter directly accreted from the background, but this 
fraction varies widely between individual haloes. Moreover, 
these fractions are little dependent on n and r^hreshi ^ 
be seen in the statistics listed in Table ^. 

2.2.2 Comparison with Other Methods 

What aHvantagps anH rlisaHvantagps Hnps this mptVinH havp 



Cole (1993)] is clearly seen. 

The absolute normalisation of the PS functions is only 
about 1 — 3 times that of the r^jjj.ggjj = 5 mass functions. 
Since the PS formula is not intended to model r^]^^^^^ — 5 
or ''thresh ~ 1000 haloes, better agreement should not be 
expected. 

Indeed, Fig. ^ shows that the haloes detected at 
''thresh ~ 1000 have a considerably lower norm alisation than 
that of the PS formula, but similar to that of Bond & My- 



relative lto others, e.g.. the BBKS peaks formalism method H 's prediction (rescaled to galaxy scales) for ellipsoidal in- 



used in ILacey et al (1993) or the peak-patch formalism of 
Bond fc Myers] (1996)? 



The advantages are that the simplifications made in the 
semi-analytical methods are not made in the N-body simu- 
lation. For example: 

(i) Most of the semi-analyti cal methods assu me spherically 
symmetric collapse (though Bond fc Myers| 1996 also con- 
sider the collapse of homogeneous ellipsoids). Warren et 
|aL| (1992) show that the haloes in the simulations analysed 
in this paper are in fact triaxial^ 



(ii) Only the most recent (e.g., Rodrigues & Thomas 1996) 



methods include spatial two-point auto-correlation function 
information. Again, this is automatically included in the N- 
body derived merging history trees (see ^2.2.4 for halo cor- 
relation functions) . The same applies for higher order corre- 
lation information. 

(iii) Characteristic dynamics of halo mergers such as tidal 
tails and particles thrown out into low-density "atmo- 
spheres" (§2.2.1) are modelled in the N-body simulation but 



not taken into account in semi-analytical models. 

Disadvantages include scientific problems such as reso- 
lution limits and practical problems such as managing large 
data files which occupy large amounts of disk space. 

A first order illustration of the similarities and differ- 
ences between this method and others is obtained by ex- 
amining the mass function plots ( Figs. The Press- 
Schechter (PS) formula derivation of Lacey & Cole (1993), 

using their Eqs.(2.11) and (2.1)|T| provides a simple analyt- 
ical description of both the shape and the time evolution of 
the mass functions. PS functions for two early and two later 
time steps are shown for comparison in each mass function 
plot, and the values of Sco = Sc{z — 0) are increased (from 
the default Sco = (5vir = 1'686) in some plots in order to show 
ways in which the PS formula could be used to interpret the 
mass functions. 

For the n — —2 mass functions, our results bracket 
the semi-analytical model results for CDM initial fluctua- 
tion spectra, which have r oughly this slo pe on galaxy scales. 



Shown on these plots are Lacev fc Silk's (1991) ma ss func- 
tion (Fig. 1; collapsed peaks) and Bond fc Myers 's (1996, 
Fig. 12) mass functions. (The latter are calculated for clus- 
ters, so masses and as are rescaled according to a length 
scahng L L/20.) 



t The factor of 2/20 in Eq.(2.1) of [Lacey fc Colcj (1993) should 
read 3/20. 



ternal halo dynamics. Since the N-body derived haloes anal- 
ysed here have a variety of triaxialities, it is not surprising 
that this provides the be st agreement. However, as remarked 
upon by Bond & Myers, their method may not be ideal for 



small haloes, since the low mass end "must find the nooks 
and the crannies" left over from the analysis devoted to high 
mass haloes, so an N-body method with sufficient resolution 
may be better for obtaining the low mass slope of mass func- 
tions. 

Lacey & Silk's (1991) peaks formalism based mass func- 
tions are somewhat steeper than the present-day PS mass 
functions and do not provide the optimal fit to our mass 
functions. 

Although the early epoch high-mass ends of the N-body 
mass functions for r^j^j-ggj^ = 5 (Figs ^, ^ are similar to the 
PS cutoffs (for Sco = <5vir)i the corresponding cuvres for the 
high density haloes (f^ijj.ggi^ — 1000) can only be matched to 
PS curves by increasing 5co (Figs ^ |^. Since 5co signifies a 
critical overdensity (according to a linear growth rate), it is 
not surprising that the high density halo high-mass cutoffs 
can be fit. (However, this would not enable the PS formulae 
to fit the later time steps, since the normalisation increases 
in proportion to Sco, which would worsen the disagreement 
for the high density haloes relative to the PS prediction.) 

Of course, what is potentially most interesting for ex- 
plaining the flatness in observed ( high surface brig htness 



galaxy) luminosity functions (e.g., Efstathiou et al 



is the agreement in the N-body derived (r^-j^j-ggj^ — 5) 
and PS slopes for n — 0. (The PS function has the form 
logj^Q dn/dM = [(n — 3)/6] logj^Q A/+const, for haloes of mass 
M). An initial fluctuation spectrum slope of ti = implies 
a less steep low mass (faint) end of the mass (luminosity) 
functions than for n — —2. 

More subtle merging p roperties, which can be compared 



with [Kauffmann fc Whitej 's (1993) PS derived Monte Carlo 
simulations, are the mass ratios of haloes and their most 
massive progenitors. The statistics of these ratios can be 
used to constrain the relation between galaxies and haloes 
via implications for the survival of disk galaxies. 

Fig. y shows that the ma ss ratios of merging haloes 
found in |Kauffmann fc White 's (1993) PS derived Monte 
Carlo simulations are similar to those found in our analysis 
for haloes detected at r^-j^j.ggj^ = 5, as for the mass func- 
tions. Since the conditions of the simulations are somewhat 
different (PS vs N-body, CDM vs n = "2), the agreement 
suggests that these ratios are quite robust with respect to 
different modelling techniques. 

Comparison with Fig. Q shows that at any given time, 
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Figure 6. Mass ratios of halo at final time stage (Mt) and the 
masses of the first and second most massive progenitor haloes 
(Ml, M2 rcsp.) of that halo at a redshift 2, for n = —2 and 
''thresh ~ ^' Circles (crosses) are for haloes more (less) mas- 
sive than 10^^ Mq. The dashed lines outline the region in the 
diagram covered by the semi-analytical haloes of Kauffmann & 
White[(]l993) for a CDM initial fluctuation spectrum and a bias 
factor of 6 = 2-5 [logiQ(Aft/Mi) values are li nearlv interoo 



lated/extrapolated in z 
White|~| 



from those in Fig. 4 of 



Kauffmann & 



the fraction of a high density halo which has already col- 
lapsed or accreted is statistically lower than for low density 
haloes. 



Figure 7. Progenitor mass ratios as per Fig. ^ for J'thresh 
1000. 



2.2.3 Merging History Trees 

Merging history tree plots (Figs ^-|ll|) are obtained by choos- 
ing a range of haloes at the final time stage and tracing back 
all the progenitors of these haloes. The line segments joining 
the circles are the key feature of the plots. These indicate 
that the halo at the earlier time is considered to merge into 
(or be "identical" to) the halo a t the l ater time according to 
the 50% criterion explained in ^ 2.1.3. The aim of the plots 
is to show connectivity over time. So, the horizontal axis 
is designed to separate the haloes according to their future 
merging activity. It doesn't directly indicate space positions, 
although there should be some correlation between how close 
two haloes are in the plot and how close they are in space 
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Figure 8. Merging History: n = —2, »"thresh ~ ^' haloes 1 — 5. 
This and the following plots show haloes detected at different 
points in space- time connected according to the criterion de- 



Figure 10. Merging History: n = —2, ''thresh ~ ^' haloes 50—60 



scribed in §2.1.3, i.e., showing which haloes merge into which. The 
horizontal axis separates individual haloes, while the vertical axis 
indicates time/redshift. (Negative redshifts indicate future times.) 
Circles indicate haloes, with radii a logarithmic transformation of 
the haloes masses (for display purposes, the specific transforma- 
tion differs between separate plots) and line segments indicating 
the merging connections. The haloes at the latest time stage, and 
the set of predecessors of any halo, are ordered by mass decreas- 
ing to the right. Numbering on the horizontal axis indicates the 
maximum number of haloes in the figure for any time stage. 
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Figure 9. Merging History: n = —2, f^jji-gsh = 1000, haloes 
1-5 



(since haloes need to be close in order to merge together 
later on). 

Much information on the merging process is represented 
in these tree plots. However, they do not show the entire 
complexity of the merging process (or "graph" in mathe- 
matical terminology). Since the halo merger history trees 
presented here start with a range of final haloes and trace 
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Merging History: n = —2, i"t}jresh ~ ^' haloes 400- 



backwards, haloes which have no descendant at the final 
time stage are not shown. As summarised in Table a sig- 
nificant fraction of the haloes at a given time stage have 
no descendants at the following time stage. A simple ex- 
planation is that a large fraction of a halo can evaporate 
in the merging process and contribute less than 50% of the 
mass of the final halo, in which case the merging/identity 
criterion adopted fails to find a descendant. Typically, about 
25% ± 25% of a halo evaporates (e.g., in tails) or forms low- 
density "atmospheres" in N-body simulations of A*' ~ 2 in- 
teracting galaxies (Quinn, 1992). These tails or atmospheres 
are likely to fall below the density detection threshold. 

Figs ^ and ^ show that merging ratios of up to around 
10:1 occur for many of the most massive haloes at low red- 
shifts, while as indicated by the maximum number of original 
haloes for any final halo, the merging ratios between early 
time stages can be much higher, as high as a few hundred 
to one in several cases for r^j^j-ggj^ — 5. 
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The difference in detection thresliolds appears clearly in 
the difference between the early time stages of Figs ^ and ^ 
At the nominal redshift z = 11, many low mass haloes have 
reached the turnaround density and are thus detected above 
'"thresh ~ 5, after which they merge rapidly. In contrast, 
haloes detected above fj-j^j-ggj^ = 1000 are mostly detected 
at z — 1-5; the few which form earlier do not merge in ratios 
as high as for the r^j-u-ggj^ — 5 haloes. Depending on the 
assumptions one wishes to make regarding the relationship 
between galaxies and haloes, either of these limiting cases 
could be interesting. 

For merging history trees of the haloes having lower 
masses at the final time stage (Figs ^ and ^l|) , very little 
merging occurs apart from the earliest few time stages. In- 
deed, Fig. |ll| shows that many of the smallest haloes either 
have only recently collapsed or are unmerged objects which 
have formed well after the first time stage. This can be gen- 
eralised by stating that the larger a galaxy halo is, the more 
original haloes it is likely to have been created from, and at 
any time in general, the more massive a galaxy halo is, the 
more haloes are likely to be merging into it. (This effect can 
also be seen in Figs g and |^.) 

That the lowest mass peaks may form fairly recently 
or may form early yet undergo no merging at all is to 
be expected in a "hierarchical galaxy formation" scenario. 
Bottom-up gravitational collapse does not only imply that 
small haloes form first and successively merge to form more 
and more massive haloes, but also that some low mass haloes 
(e.g., from low amplitude short length-scale perturbations, 
which must exist if the initial perturbation amplitudes are 
part of a zero-mean Gaussian distribution) continue to form 
at late epochs. 

The recently forming small haloes in our halo merger 
trees could be used to model the (observed) existence of 
low mass, young, low redshift galaxies, such as the dwarf 
Spheroidals. An interesting case is the SBS 0335-052(W) 
pair of dwarf galaxies at a redshift of z = 0-014, which have 
a (stellar) age of not more than 10* yr an d have a met allicity 
(0/H) around 1/40 of the solar value ( [[zotov et al] 1997). 
These two dwarfs are extremely hard to understand as any- 
thing other than young, low redshift, low mass galaxies. For 
such galaxies to provide a constraint on galaxy formation 
models, precise estimations of number densities and age dis- 
tributions would be useful — though obviously difficult to ob- 
tain without strong observational selection effects. 
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Figure 12. Spatial two-point autocorrelation functions of haloes 
in n = —2, ''thpesh ~ ^ model, shown as logiQ[^(r)] against 
logjQ(r), where r is the comoving halo pair separation in Mpc. A 
solid line of slope —7 = —1-8 is shown for comparison. 




Figure 13. Evolution of the spatial two-point autocorrelation 
function, = C(5 h~^Mpc, z) vs redshift (1 -|- z). (The highest 
redshift correlation function for n = 0, '"thj-esh ~ 1000 is too 
noisy to deduce a ^0 value.) Also shown are lines fitted to all but 
the highest z point for n = and to the four points with lowest 
z for n = —2. 



2.2.4 Halo Correlation Functions 

Another important statistic of the haloes is their spatial two- 
point auto-correlation function, (,{r). The natural inclusion 
of ^ in N-body simulations is something not usually p resent 
in semi-analytical galaxy formation modelling. Indeed, Vane 



et al 



s (1996) application of Jedamzik's (1995) approach to 



the " cloud-in-cloud" prob lem of the Press- Schechter formal- 



ism ( Press fc Schechter 1974) shows that the mass func- 



tion of the l atter suffers significant ly from the lack of in- 



clusion of ^. [Nagashima fc Gouda (1997) confirm th 



semi-analytical me rging history tree simulations using Ro- 



drigues fc Thomas 's (1996) modification of the Block model 
( |Cole fc Kaiseil 1988). 

Fig. I12I shows the correlation functions for the haloes de- 



tected for n — —2, r^-[jj.ggjj = 5. These functions have slopes 
—7 ~ —1-8, where 

-{3-|-e-7) 



^^{ro/rni + zy 

r and ro are expressed in comoving 
parametrises the growth rate of ^ o(z) 



either from observation or theory (Groth & Peebles 1977) 



(1) 

coordinates, and e 
= £,(ro, z) , deduced 



The slopes of ^ for the other combinations of n and 
''thresh ''■''^ similar, so are consis tent wit h tha t observed for 
' <1_8 (e 
1992) 



galax ies, i. e., 1-7^7 



1983; Loveday et al 



Peebles 1980 



Davis & Pee- 



bles 

that 7 evolves with time ( Infante 
ers do not detect this 



Some observations indicate 
1995), but oth- 
( iHudon fc Lilly! 1996)- N-body sim- 



ulation derived correlation function slopes vary around the 
same values, but are sensitive to whether the correlation 
is that of particles (identified as galaxies), of haloes or of 
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Table 7. Low redshift spatial correlation function growth rates 
e (see Eq. |^. Only the four lowest redshift points are used for 
n = -2. 



'"thresh 


n = 


ra = -2 


5 


-0-6 ±01 


-0-6 ±0-6 


1000 


-0-4 ±0-3 


-0-3 ±0-3 



mass-d ensity ( Davis et al. 1985; [Efstathiou et al 



Sugiiiohara 1991; suginohara & Sutc 1991; Melott 1992) 



1988; 3utc 



The haloes of our N-body derived merging history trees 
therefore have reahstic correlation slopes and do not suffer 
from the lack of a correlation function of the traditional 
semi-analytical methods. 

Another interesting property of the correlation function 
is the growth rate of ^o, parametrised by e. This is shown in 
Fig. h_3. Since the simulations used here are not normalised 
conventionally, the values of ^o{z ~ 0)^ are less useful than 
the change of with redshift. 

To the extent that Eq. ^ is a good approximation, direct 
observational estimates of e at varying ^ mpH (median red- 



shift) include e = 1-6 ±0-5 (Warren et al 



( Shepherd et al. 



Le Fevre et al. 



1993, 



e = -2-0 ± 2-7 ( pole et alj 1994, z^^^= 0-16) 



^med 



1997, 
(1996) 



0-4), 
e ^ 2-2 

2ined= 0-37) and e = 2-8 [cf. §4.1(3), 
adopting ro = Mp c and 

1992, 



cZj^gj(Stromlo-APM)= 15, 200kms ^(Loveday et al. 
1995) 

1 ^medlCFRS)— 0-56]. By hypothesising major changes 
in visible galaxy populations from low to high redshifts, 
these values of e could be lowered. 

Simple theoretical values of e include those expected for 
clustering fixed in comoving coordinates (e = 7 — 3~— 1-2, 
so that the index in Eq. ^ is zero) ; clustering fixed in proper 
coordinates on small scales ["stable clustering", e = 0; since 
the numbers of "clusters" changes by (l-|-z)^ and the number 
of galaxy pairs by (1 -I- z)'' , the factor in Eq. ^ for r in 
proper coordinates is (1 -I- z)~^]; and linear growth of initial 
fluctuations for an f^o = 1 universe (e = 7 — 1 ~ 0-8) . 



More detailed a nalyses include those of Peacock (1996) 
and Matarrese et al. (1997), who discuss departures from the 
power law of Eq. |l| bias factors and difi'erences between the 
linear, quasi-linear and non-linear overdensity regimes. In 
particular, the transition between the linear and non-linear 
regimes can give r ise to e ~ 2-8 [for "/ — 1-8 in Eq. (29) of 
Peacock fc Dodds| (1996) and f^o = !]• 

By contrast to the above estimates, the values of e for 
our haloes are those for which haloes merge between time 
steps. This means that as haloes approach each other and 
merge, halo pairs are replaced by single haloes, so that the 
correlation decreases more slowly than if the haloes' iden- 



tities wieie kept sepai 
derived^ here (Table W 



ate. Because uf this, the values uf t 
\ lie between those expected for clua — 
tering flxed in comoving and proper coordinates respectively, 
and are lower than the more precise of the observational es- 
timates. 

Implications of the effect of merging, in particular on 

t These could be reinterpreted via a change of time or length 
scale (earlier time or larger length scale), or the region studied 
could be considered to be a denser than average region in a low 
density (e.g., Co = 0-1, Aq = 0) universe model. 



=1 



12 3 4 

time to merger in Gyr 

Figure 14. Method by which time steps could be interpolated. 
Histogram of time of radial infall of haloes (at t = 9-3 Gyr for n = 
—2, '"^iij.gsjj = 1000) considered as point particles which fall into 
isothermal potentials of the combined masses (at t = 12-3 Gyr) 
of the multiple merger products. 



the angular c orrelation function, have been presented in 
more detail in Eloukema & Yoshii (1993). 

Another property to note is that for the n = — 2 model, 
^0 initially decreases until z ^ 1, after which it increases 
to the present. This is because during the transition from 
linear perturbation amplitudes to non-linear collapse, small 
length-scale perturbations on top of large length-scale per- 
turbations have their density boosted, and so can collapse 
well before small length-scale perturbations in low density 
regions. During this period, perturbations in the low-density 
regions are not represented by collapsed objects. Hence, the 
mean number density of collapsed objects used to normalise 
the correlation function represents the numbers of haloes 
in only the high-density regions, but includes the volume 
in both low and high-density regions. The mean density is 
therefore biased low, so the correlation function amplitude 
is biased high. As a larger fraction of perturbations collapse, 
this initially high bias in ^ decreases until normal correlation 
growth takes place. 

This effect is not seen for n = because large length- 
scale perturbations have relatively less power for n = 
than for n = —2. Possible consequences of s uch a "Dec reas- 
ing Correlation Period" (als o reported by Roukema 1993 
and 



Brainerd & Villumsen 



1994) on observati ons of the 
angular correlation function are discussed in 3gawa et 
ai] (1997). From Fig. ^ the rate of correlation decrease 
is d(logio Co)/d[logio(l + z)] ~ 0-3 ± 0-15 to z ~ 0-5 - 1. 



2.2.5 Merging Between Output Time Steps 

The merging history trees of haloes could in principle be 
analysed with a higher time resolution than in the figures 
presented here, merely by using many more output time 
steps from the N-body simulation. The only drawback would 
be practical handling of computer disk space. A single time 
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step for the Warren et al. (1992) simulations normally oc- 



cupies 32 Megabytes. 

Here we note an alternative to using more output times. 
Interpolation between time steps, making analytical as- 
sumptions for the trajectories of the haloes, can give in- 
terpolated merger times which are statistically consistent 
with the more accurate trajectories calculated in the full 
N-body simulation. Fig. |l4| shows that the approximation 
of point particles falling into isothermal potentials gives an 
approximately decreasing merger rate, which is consistent 
with the overall merging rate. A small fraction of objects 
have estimated merger times greater than the output time 
step interval, i.e., greater than the time required for a merger 
according to the full N-body simulation. This is a limitation 
to (at least this) analytical interpolation, but the fraction of 
such cases is small enough that this interpolation technique 
may be a fair approximation statistically. 

However, N-bo dy an alyses of few-galaxy sys tems (e.g., 
Prugniel & Combet 1992; Hernquist & Weinberg 1989) sug- 



gest that merger times may be difficult to predict by any 
simple analytical method, so caution is needed. Although 
the results should not differ significantly from those of the 
authors just cited, use of such interpolations for the calcu- 
lation of N-body merger history trees could be verified by 
comparison with the intermediate time steps of the N-body 
systems. Interpolation is not adopted for this paper. 



3 GALAXY FORMATION APPLICATION 
(STAR FORMATION AND EVOLUTION) 

The example of a galaxy formation recipe adopted here com- 
bines (a) gravity: halo merging history trees are derived from 
density peaks detected above a given overdensity thresh- 
old in 1-1 million particle pure gravity N-body simulations; 
and (b) other physical processes: one galaxy is assumed to 
form in each halo; we insert an observationally inspired star 
formation rate due to merger-induced starbursts (and/or a 
quiescent star formation rate); and we combine this with 
the initial mass function and stellar evolution modelled in 
Bruzual's (1983) evolutionary population synthesis code. 



3.1 Modelling Starbursts to Occur on Merging 

Star formation physics is represented by an initial mass func- 
tion (IMF) and a star formation rate (SFR). The separabil- 
ity of these two functions appears to be a reasonable as- 
sumption according to both observation and theory, e.g., 
the recent ISO (Infrared Satellite Observatory) observations 
of both infrared-ultra-luminous and "normal" star-forming 
galaxies are consi stent with exis tence of a uni versal IMF 
( Lutz et al.| 1996; [Vigroux et al] 1996; |Gilmore| 1996). The 
IMF is discussed in ^3.2| , while this section describes the 
modelling of the SFR, based on starbursts. 

For the purpose of this first-order examination of the ef- 
fect of merger-induced starbursts on galaxy luminosity evo- 
lution, we only use very simple models of the starbursts. 
Observational and theoretical motivations are as follows. 

An early observational model of a starburst (not neces- 
sarily caused by a merger), is that of Rieke etal. (1980). 

Rieke et al. (1980) model the starburst in the nucleus of 
M82 via evolutionary population synthesis. They find that 



instantaneous (Dirac delta function) and constant star for- 
mation rate models fail to produce the observed spectrum, 
but exponential decay SFR models with an IMF with a low- 
mass cutoff well above a solar mass are necessary. Their best 
models (say, D and F) have the e-folding time in the SFR 
to = 2xl0'^yr and to = IxlO^yr and run for t = SxlO^yr and 
t — l-6xl0^yr respectively. Both have IMF's with a = 2 and 
the mass range 3-5 — 31M0. The mass turned into stars is 
~ 1-5 — 2xlO*M0, this being constrained to be less than the 
total mass in the nucleus, estimated as 3x10*Mq by Rieke 
et al. This constraint is the major reason for the need of a 
high low-mass cutoff. If a normal low-mass cutoff is used, 
then the mass actually present in the nucleus is insufficient 
to generate the observed luminosity. 

Rieke (1991) describes more recent observational con- 
straints on the models for M82, finding that the above con- 
clusion still holds using more modern stellar evolutionary 
tracks in the models. 

Scoville & Soifer (1991) argue from IRAS far-infrared 
data that "virtually all of the strong global starbursts occur 
in ... starburst-infrared galaxies," where the latter are de- 
fined as "those with Lm/Mn^ significantly higher than in 
normal galaxies," and that starburst-infrared galaxies corre- 
late highly with "the occurrence of a recent [galaxy-galaxy] 
interaction." They argue that such global starbursts require 
the progenitor galaxies to be of comparable mass in order 
to generate such activity. 

While this result doesn't necessarily imply the converse, 
i.e., that all mergers of similarly sized galaxies induce major 
global starbursts, the converse is a fair hypothesis. With 
the assumption of this converse, Scoville & Soifer find that 
the high luminosity end of the infrared (galaxy) luminosity 
function from the IRAS survey is consistent with 0-2% of all 
spiral galaxies undergoing global starbursts at the present 
with lifetimes equal to the dynamical times of large galaxies, 
~ 1 — 2xl0®yr. For the most IR-luminous galaxies they find 
SFR's of 10 - lOOMoyr"^ They don't find a high low-mass 
cutoff necessary for their models to fit the observations. 

Norman (1991), citing the models of merg ing gas-rich 
disk galaxies of Hernquist and Barnes, (Barnes 1990, equal 
mass galaxies; Hernquist 1989, differing mass galaxies), de- 
scribes a qualitative three-phase model to take into account 
gas falling into the galaxy centre. The three phases essen- 
tially correspond to proportions of the gas having fallen to 
the centre. Published star formation models following these 
three phases separately are not cited in the article, but would 
obviously be of interest. 

Norman (1991) also argues that constant SFR models 
of starbursts satisfy an observed sparsity of post-starburst 
galaxies relative to starburst galaxies, but that instanta- 
neous SFR models predict too many post-starbursts. 

Hence, for this example application of an N-body de- 
rived halo merger history tree, a starburst with a constant 
SFR is chosen. 

For pairwise mergers, the following canonical values 
are chosen. We normalise the rate of the starburst for a 
"typical" large galaxy merger product to be an SFR of 
■i/'o = SOMeyr"^ as in the models of Scoville & Soifer (1991). 
The low-mass cutoff in the IMF used here (0-08Mq, see 
Eqn ^ is consistent with Scoville & Soifer's value of O-IM©. 

We consider the merger product to be the remnant of 
two large galaxies each of gas mass M^'^ = 10'°AfQ, to- 
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tal mass Mq'^'" = W'-'^Mq and halo radius 50kpc 
a dynamical time t^yn = (Gpo)"^''^ ~ 2 



This 

gives a dynamical time tj-yn --i/o 
where the mean density of either galaxy to its halo radius 
is Po = SxlO^'Mokpc"^. The modelling by |Barne^ (1990), 
Hernquist (1989) or earlier non-gaseous N-body models such 
as those of Quinn & Goodman (1986) find that the merger 
takes place over only a few dynamical times. So the choice 
of progenitor galaxies here gives a t^y^ matching Scoville & 

Soifer's burst duration of fg"'^^* ~ 2xl0*yr, which is chosen 
as the canonical burst duration. 

In this canonical case, 50% of the total gas mass is used 
up in the burst. To sum up, we have 



i>o = SOMoyr 



M, 



halo 



10^° M, 

12 i 



= 10'^ M, 



Po 
.burst 
to 



= Ox 



lO'^Afokpc"^ 



2x10 yr. 



(2) 



The canonical values are then scaled for haloes of arbi- 
trary mass. To first order it seems reasonable that the ki- 
netic energy available for generating star formation should 
be proportional to the mass of the smaller halo and that the 
SFR should also be proportional to the total amount of gas 
available. So, the SFR is scaled by the mass of the smaller 
halo times the ratio of combined gas mass to 2Mq^. 

Since we consider the duration of the starburst to be the 
order of a dynamical time, t'^'^''^* is scaled by p~^^^, where 
p is the density of the larger halo. 

Given a coarse time resolution in the merging histories, 
each merger is identified by the code as a multiple merger — 
e.g., seven haloes merge to one — instead of as a series of 
several individual pairwise mergers. If we consider the ap- 
proximation that each of the pairwise mergers involves the 
largest progenitor halo, then a compromise can be carried 
out as follows. 

(a) Have a single burst with the above normalisation. 
Scale the SFR by the sum of the masses of each of the smaller 
haloes (i.e., all but the largest) and by the ratio of the com- 
bined gas mass from all progenitor haloes (i.e., including the 
largest) to 2M^'^ , giving 



■ip{t) = i>o 



(3) 



where tp{t) is the star formation rate (in Moyr'^), the pro- 
genitor haloes are labelled by i, and i™^^ is the label of 
the progenitor halo of greatest mass. (For the group-finding 
algorithm used in this article, i™^^ ~ 

(b) Scale the starburst duration as above, by p~^^^, 
where p is the density of the largest halo, giving 



.burst .burst 

t —to 



Phalo 
Po 



-1/2 



(4) 



Two modifications may have to be applied in some situa- 
tions. Firstly, where the starburst at this rate uses up more 
gas than is actually available, it is truncated at the point 
of time when zero gas mass is left. Secondly, if the duration 
of the starburst is longer than the time interval to the next 
time step, it is truncated at that next time step. 

While this modelling of multiple merger- induced star- 



bursts with large time steps may make the luminosity evo- 
lution more temporally discrete than it should be, it does 
conserve physical quantities in line with the observational 
and theoretical constraints discussed above, and should be 
sufficient for this demonstration of the use of N-body based 
merger history trees. 

3.2 Connection with Stellar Evolutionary 
Population Synthesis 

We use stellar evolutionary population synthesis to combine 
star formation and stellar evolutionary physics. A version 
of Bruzual's stellar evolutionary population synthesis code 
(Bruzual 1983) which is essentially that of 1983, but with 
some updating and conversion to Unix, is th e primary pop- 
ulation synthesis c ode used, but the code of 



Rocca-Volmerange (1987) and Rocca-Volmerange & Guider- 



Guiderdoni & 



doni (1988) (1993 version) gives similar results. In the inter- 
actions of the programs presented above with this code, the 
SFR history and masses of galaxies and gas are determined 
outside of the Bruzual routines or by amended versions of 
the Bruzual routines. The return of gas from supernovae to 
a galaxy was turned off for test purposes but otherwise left 
on. The loss of this supernova gas from a galaxy was not 
invoked, neither was the option allowing infall of gas to a 
galaxy. 

The initial mass function (IMF) used was the default 
IMF chosen by the code, after Scalo (1986) (e.g., Fig. 16 in 
Scalo, 1986). Where /(m) oc m"'^^^^ is the number of stars 
born per unit (linear) mass in a given mass range, the IMF 



used are 




-2-60, 


0-05 < M < 0-08Mq (brown dwarfs) 


-2-60, 


0-08 <M< 0-18Mq 


0-01, 


0-18 <M < O-42M0 


1-75, 


0-42 <M< O-62M0 


1-08, 


0-62 <M< 1-18M0 


2-50, 


1-18 <M< 3-5M0 


1-63, 


3-50 <M< 75Mq. 



(5) 



The Bruzual code normally works by using simple ana- 
lytical expressions for the SFR history, so that no numerical 
effects (e.g., rounding errors) can be introduced at this stage. 
Numerical effects can of course be present when the galaxy 
spectral energy distribution is calculated, since only finitely 
many points representing different stellar ages are present 
for each of the finitely many mass tracks. 

However, this has the disadvantage that one cannot sim- 
ply stop the code after a certain time step, save the popu- 
lation data, start up the program from scratch, read in the 
saved population data and continue on as if the program 
had never stopped. The population data could be stored 
and later read back in, but this would round each star's age 
to the appropriate stellar evolution track age at every time 
step, making cumulative errors. 

Bruzual's code was therefore modified slightly in order 
to allow use of numerical SFR histories. For each time step 
and each halo, an array of time points from that time step 
to the following time step and the corresponding array of 
integrals of the SFR are stored. These integrals of the SFR 
are the total number of stars created since the first star 
formation in any of the progenitor peaks which end up in 
the present peak being worked on. Because the integral of 
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the SFR is used, the errors are not cumulative, and in the 
present version are ~ 0-1%. 

The program which applies stellar evolutionary popu- 
lation synthesis to the merging histories stores an SFR his- 
tory for each peak as it is evolved to the next time step, 
adds these together for merged peaks, and from that point 
on evolutionary population synthesis applies just as it would 
for an isolated galaxy, except that the merging history tree 
information is contained implicitly in the complex shape of 
the SFR history. 

For the present model, population synthesis is applied 
by optionally having an exponentially decreasing SFR be- 
tween mergers and optionally having starburst SFR's com- 
mencing at each merger. Gas masses and total masses are by 
default conserved, i.e., the gas and total masses of a galaxy 
are the respective sums of the gas and total masses of pre- 
decessor galaxies, except that matter accreting directly onto 
the density peaks is added as gaseous mass. If both exponen- 
tial and burst star formation are turned on simultaneously 
(probably the most realistic model) the SFR's are simply 
added together, conserving the number of stars created. 

3.3 Luminosity Functions 

This combination of N-body generated merging history trees 
with other elements of a simple galaxy formation recipe is 
sufficient to generate luminosity functions which are compa- 
rable to the estimated present-day luminosity function, and 
so could be extended to an exploration of galaxy formation 
parameter space, in which the parameters are free rather 
than motivated by observation as in the present work. 

The merging history trees presented above are derived 
from power law initial spectrum N-body simulations in- 
tended for comparison of the relative properties of haloes of 
different masses, so need to be renormalised in a cosmologi- 
cal context. For simplicity, an increase in the length scale by 
a factor of two is the renormalisation adopted. Masses scale 
as the cube of the length in order to leave the numerical op- 
erations in the N-body simulations and the merger history 
tree analysis unchanged and time is not rescaled. This brings 
the number density of the haloes close to that of observed 
galaxies. 

A baryo n fraction of 1 0% is adopted. Nucleosynt hesis 
results (e .g 



t (yr) 



Figure 15. Star formation rate (SFR) history for tlic most mas- 
sive galaxy at the final time step in the n = —2, ^'thresh ~ ^' 
exponential-f burst model, fj. = 0-15. The SFR in stars yr"'^ is 
plotted against time since the formation of the first progenitor of 
the galaxy. 
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1995; [Kundic et al.| 1997) make this a reasonable Figure 16. A^iumZ-^IIIaJ ^^^^ ^ame values of ^niaj 



round number for a critical density universe with a null cos- 
mological constant. [Some other estimates of similar baryon 
fractions include Galaxy estimates of 0-7 out to a Hol mberg 



radius, or down to possibly 0-07 for the whole Ga laxy (Free 
1987); 0-25 for galaxy clusters (e.g., | arazin] 1987); and 



OT over the whole range of cos mological str uctures from 
dwarf spheroidals to rich clusters (Blumenthal 1988).] 

All baryonic matter is assumed to be potentially star- 
forming material. 

An example of an SFR history (for the galaxy resident 
in the most massive final halo in the n — — 2, r+hresh = 5 
exponential-plus-burst model) is shown in Fig. Il3. Both the 
exponentially declining "quiescent" star formation rates and 
the merger-induced bursts are clearly visible. Bruzual's SFR 
parameter has the value = 0T5 [fi is the proportion of gas 
turned into stars in a (non-merging) galaxy within 1 Gyr] . 

Mass-to-light ratios, A^ium/^IIIaJ' (using rest frame 



a run on the n = — 2, '"thresh ~ ^ model merging history with 
exponential and burst evolution turned on and Bruzual's SFR 
parameter fi = 0-15. The masses detected at the time stage ti 
and the luminosities resulting at the end of the interval [ti,ti^i) 
are used to obtain A^lum/^niaJ values labelled ti in this figure. 



values of Lm^j, i.e., no K-corrections) for the galaxies in 
this model at different time steps are shown in Fig. |l6[ 
These mass-to-light ratios are somewhat high relative to op- 
tical galaxy estimates (e.g., 2 < A^disk/^s^isi^^ ~ -^0-^0^ 
derived from the inner part of the optical/H I rotation 
curves of disk galaxies without bulges, Freeman 1987; stel- 
lar population values, Larson 



Tinsleyl 1978) but simi- 



lar to those estimated from X-ray emission in ellipticals 



(20 < Mtnf /L B < 30 Mr.-,Lr^ within a radius r 



40 kpc, Canizares 1987) 
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Figure 17. IllaJ luminosity functions for the n = —2, ^^j^i-ggj^ = 
5 model, exponentially decaying plus burst SFR, Bruzual's SFR 
index fi = 0-15. The model curves are of the line styles indicated, 
labelled by the beginning of each time stage interval, and com- 
pared with a Schechter function with locally estimated param- 
eters (solid line, Eq. Luminosities are expressed in absolute 
magnit udes. A^jjiaJ' observer's frame, i.e., K-corrected 



(Wirtz 1918); densities are in logj^Q(N Mpc ^ mag ^) in comov- 



ing coordinates. The following plots show the same quantities 
for different n and '"thresh- 
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Figure 18. IllaJ luminosity functions for n = —2, ''thresh ~ ^ 
model, (exponential decay only SFR), Bruzual's fi = 015. 



The corresponding luminosity functions for the galax- 
ies in this model are shown in Fig. while luminosity 
functions for either exponentially decaying or burst SFR's 
only (not both) are shown in Figs |l^ and |l^ (for n = 
—2, ''thresh = 5). A local luminosity function parametrised 
as a Schechter function 



dN/dL = <;/)*(i/L*)" exp(-L/L*)d(i/L*) 



(6) 



( Pchechteij 1976) where 0* = l- 56xl0~^/i^ Mpc" '^, M* = 
— 19-6 -I- Slog jg /i and a = —1-1 (Efstathiou et al, 1988), is 
shown for comparison. [This estimate is similar to the more 
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Figure 19. IllaJ luminosity functions for n = —2, ''thresh ~ ^ 
model, burst-only SFR. In this and the following plots of the lu- 
minosity function for burst-only models, the luminosity functions 
at the first time stage are missing. This is a property of the model: 
stars are only formed when mergers occur; this first occurs at the 
second time stage. 



the CfA2 estimate of Marzke et al 



recent Stromlo-APM estimate of Loveday et al. (1992), while 



(1994) differs from both 
in having a value of M* about 0-7 mag fainter.] 

The luminosity functions for the n = — 2, ^thresh ~ 5 
exponential-plus-burst and exponential-only SFR's are sim- 
ilar to one another (Figs |l8|) , as the rate of the exponen- 
tial SFR alone is enough to use up most of the gas, leaving 
little possibility for the bursts to make any difference be- 
tween galaxies of differing merging histories. The luminosity 
function slopes are steeper than the observational faint end 
slopes (for field galaxies) , and are similar to the steep slopes 
expected from the mass functions (Figs ^ - §). The other 
three combinations of n and '"thresh ^r the exponential- 
plus-burst and exponential-only models give similarly steep 
luminosity functions. This re sult is similar to that of other 
semi-analytical models (e.g., [Kauffmann et al" 1993). 

The "bump" at the bright end of some of the later lumi- 
nosity functions (particularly in Fig. |l^) combined with the 
steep slope is suggestive of luminosity functions in clusters, 
but a cluster environment is one in which merger- induced 
star bursts could be expected to be the more important 
rather than the less important star formation mode, so it 
is not clear whether or not this is a useful characteristic of 
the models. 

The faint end sudden drop (at M ^ — 16 in these two 
figures) is simply due to the resolution limit of the models. 

The n — 0, ^thresh ~ 5 burst-only model (Fig. po| is 
quite striking in having a "faint end" slope almost precisely 
as shallow as the slope of the observational luminosity func- 
tion. In addition, the bright ends of the earlier of these lu- 
minosity functions have steep drops as in the observational 
luminosity function, though this occurs for galaxies a few 
magnitudes too bright. 

How significant are the luminosity functions of the 
n — 0, ''thresh ~ ^ burst-only model? For a fiat rotation 
curve for the Galaxy with density falling off as r~^, the 
detection threshold of ''thresh ~ 5 corresponds to a Galaxy 
halo radius of about 1500 kpc, larger than any claimed value. 
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Figure 20. IllaJ luminosity functions for n = 0, ^'^jjj.gg]-! = 5 
model, burst-only SFR. 
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Figure 21. IllaJ luminosity functions for n 
1000 model, burst-only SFR. 



-2, r 



thresh 



Apart from the recent observational suggestion of Honma 
& Sofu| (1996) that the halo is only 15 kpc in radius, a 
conventional estimate of rhaio ~ 100 — 300 kpc would still 
imply that the detection threshold used here overestimates 
the mass of the Galaxy by an order of magnitude [since 
p oc M(< r) oc r]. Also, members of the Local Group 

would missed, since one galaxy per halo and maximal merg- 
ing are assumed. 

These problems could conceivably be corrected by a re- 
duction in galaxy masses, which would reduce the luminosi- 
ties, and by an increase in the normalisation to account for 
the missing Local Group galaxies, whose analogues would 
presumably be detected (in the model) in isolation in other 
regions of space. However, while the former would improve 
the fit, the latter would worsen it. 

The n = 0, f^-j^j.^gj^ = 1000 burst-only model does in- 
clude these smaller galaxies ("halo substructure") and re- 
sults in a luminosity function which is slightly too steep 
(Fig. H). The n = -2 models (Figs |l|, |2^) also have slopes 
which are relatively shallow, but none as shallow as that for 
the n — Q, r^-j^j-ggj^ = 5 model. 

The reason why the burst-only models have shallower 
slopes can be explained fairly simply: smaller mass haloes 
typically undergo less merging, so become relatively underlu- 
minous relative to higher mass haloes. Star formation dom- 
inated by merger-mduced bursts could therefore provide a 
simple way of reducing the faint-end slopes expected from 
the Press-Schechter function. 

As might be expected from the similarity in the mass 
functions (Figs ^, ^) for three of the burst-only mod- 
els, the luminosity functions are also similar to one an- 
other. However, a characteristic of the mass function of the 
n — —2, r^]-|j.ggjj = 5 model appears (marginally) to be 
shared by the corresponding luminosity function: the slope 
seems to become shallower with time. If this were significant, 
it could be of considerable inter est in reconci ling faint galaxy 
counts (e.g., Tyson & Seitzer 1988; Tyson 1988) with the 
flatness of the locally estimated luminosity function, though 
it would also need to show evolution in the bright end of the 




-30 -15 
IllaJ lum of peak 



Figure 22. IllaJ luminosity functions for ra = 0, ^'thresh ~ 1000 
model, burst-only SFR. 



l uminosity functi on as dedu ced from recent redshi ft surveys 



( |CFRS-VI| 1995, 



Ellis et al 



1996 and Cowie et al 



1996). 



The evolution in slope could be attributed to the evo- 
lution in the Press-Schechter mass function shown in Fig. ^ 
with the addition of a systematic reduction in slope induced 
by our burst-only SFR formalism. 

Both the flatness of the n — 0, r^-j-^-ggj^ — 5 burst- 
only luminosity functions and the decreasing slope of the 
n — —2, f^hresh ~ ^ model suggest that these parameter 
combinations should be interesting for future explorations 
of galaxy formation models using N-body derived merger 
history trees. 



4 CONCLUSION 

The algorithms and results discussed above show that this 
method of deriving halo merging history trees directly from 
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N-body simulations is feasible and easily combined with evo- 
lutionary stellar population synthesis for exploration of a 
simple g alaxy formation model. The method has been ap- 



servations such as that of a young, l ow redshift, low metal- 
licity dwarf galaxy pair (Izotov et al, 1997) can be explained 



plied to Warren et al.'s (1992) N-body simulations with 
n = and n = —2 power law initial perturbation power 
spectra, using r^j^j-ggj^ = 5 and r^j^j-ggj^ = 1000 overdensity 
detection thresholds to detect dark matter haloes. 

Several subsets of the merging history trees have been 
plotted, directly showing quantitative patterns of halo merg- 
ing from fully non-linear calculations. The merging history 
trees were combined with simple, observationally normalised 
proportionality assumptions for star formation rates and 
stellar evolutionary population synthesis in order to demon- 
strate a simple galaxy formation "recipe" , in which star for- 
mation is either exponentially decreasing, induced by merg- 
ers, or both. In addition, interesting properties of halo for- 
mation and evolution were noted. 

The galaxy formation recipe adopted shows that if most 
star formation occurs as merger-induced bursts, then some 
flattening of the faint end of the galaxy luminosity function 
relative to that expected from the mass functions may be ob- 
tained. Indeed, for the analysis using (a) a white noise per- 
turbation spectrum slope (n — 0) and (b) a detection thresh- 
old typical of that for a perturbation just reaching the turn- 
around density (''thresh ~ 5)' the resulting present-day lu- 
minosity function has a faint-end slope similar to a = —IT, 
i.e., that estimated for local, field, high surface brightness 
galaxies. Moreover, for the same value of ^'thresh ^^i^ a per- 
turbation slope similar to that of CDM, there marginally 
appears to be a steepening of the faint end slope for increas- 
ing redshift. 

However, since this is only a simple test application of 
merging history trees to galaxy formation modelling, and 
the N-body simulations are designed for the study of relative 
rather than absolute halo properties, these results should be 
taken as indicative of worthwhile galaxy formation "recipes" 
to explore further, rather than as a definite explanation for 
a luminosity function with a flat faint end. 

Other interesting by-products from this method of 
analysing non-linear gravitational collapse are properties of 
hierarchical halo formation. 

(1) Individual merger rates can be very different from 
average merger rates and the fraction of mass coming 
from accretion can be quite high. For example, for n = 
0, ''thresh ~ 5, the mean number of haloes which collapse at 
any time stage and end up in a halo at the flnal time stage 
is 7-4, the standard deviation in this quantity is 20-7 (Ta- 
ble ^) and the maximum is 233. For the same halo detection 
parameters, a mean fraction of 32% of the mass in the flnal 
haloes comes directly from accretion of uncoUapsed matter 
(Table m, but this fraction varies widely between different 
haloes. (For the other three parameter choices, the merger 
ratios are about a factor of two lower, while the accretion 
percentages are about the same.) 

(2) Low mass haloes may either form at very recent cos- 
mological epochs or may undergo no mergers at all through- 
out a Hubble time. The former may be counter-intuitive, 
but can be simply understood as being due to the existence 
of small amplitude, small length-scale perturbations at early 
epochs, expected for a zero-mean Gaussian amplitude distri- 
bution. Some (if not many) perturbations have small initial 
amplitudes: these only become non-linear recently. Thus, ob- 



naturally within a hierarchical halo formation model. 

(3) If there is sufficiently more power on large than small 
scales, e.g., n — —2, then the first fluctuations to collapse 
(forming haloes) are only those inside of large length-scale 
perturbations, so that the spatial distribution of the haloes 
is initially much less uniform than that of the linear (and 
non-linear) fluctuations in general. In other words, the am- 
plitude of the spatial correlation function of the haloes may 
start from a highly biased value (relative to the linear per- 
turbation theory expectation), and decrease in bias until a 
transition redshift when "normal" correlation growth com- 
mences. Observational consequences of such a "D ecreasing 
Correlation Period" are discussed in Ogawa et al, (1997). 

Many further developments of the use of N-body derived 
merger history trees are obvious: different assumptions for 
processes (2)-(5) cited in the introduction could be adopted; 
a mechanism (e.g., dissipation) could be chosen to change 
(6), so that galaxies would not merge every time their par- 
ent haloes merged; or at the expense of using more computer 
disk space, N-body simulations with higher A'' or larger num- 
bers of output time steps could b e used. Use of more p hysi- 
cally complex group-flnders (e.g., Gelb & Bertschinger 1994) 



should also make some differences to detailed properties of 
the merger histories. Syste matic compa risons between recent 
analytical predic tions (e.g., Yano et alj 1 996), Merging-Cell- 
Models models ( |^odrigues fc Thomas 1996) and N-body 



derived merger history trees could be made. An example of 
a different assumption for processes (l)-(5) would be to re- 
place the observationally inspired parameters adopted here 
by free parameters which would be varied to see how sensi- 
tive the reduction in the slope of the faint end of the lumi- 
nosity function is to these parameters. 
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